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The transition metal-oxygen bond appears very prominently throughout chemistry and solid- 
state physics. Many materials, from biomolecules to ferroelectrics to the components of supernova 
remnants contain this bond in some form. Many of these materials' properties depend strongly on 
fine details of the TM-O bond, which makes accurate calculations of their properties very challenging. 
Here we report on highly accurate first principles calculations of the properties of TM monoxide 
molecules within fixed-node Diffusion Monte Carlo and Reptation Monte Carlo. 



I. INTRODUCTION 

Transition metal chemistry is an exciting area of re- 
search that has implications in fields from biological 
physics to astrophysics. Transition metals can form many 
types of bonds and transition metal solids exhibit useful 
properties like ferroelectric and ferromagnetic ordering. 
This interesting physics comes from the <i-shell states 
which are very close in energy to the outer s— states, 
along with strongly correlated electrons that make accu- 
rate first principles calculations particularly challenging. 
Benchmark calculations are particularly useful to deter- 
mine what level of accuracy one can obtain from a given 
method, although the precise bonding pattern can vary 
from system to system. 

Many authors (most recently, Refs 0,0) have stud- 
ied the transition metal monoxides using Density Func- 
tional Theory and post-Hartree-Fock methods. The per- 
formance of these methods is less reliable whenever tran- 
sition metals are included in a system. In particular, 
the calculation of dipole moments is challenging because 
it is rather sensitive to the details of calculations and 
sizes of employed basis sets. The results can sometimes 
be far from experiment. The TM-0 bond is the driving 
force behind many perovskite and earth materials, which 
have been noted as having significant errors in the unit 
cell volume within DFT[30|, while being too large for 
post Hartree-Fock methods to be applied. We have had 
some success applying ground-state quantum Monte 
Carlo (QMC) to the binding energies of the TiO and 
MnO molecules, which hinted that QMC may be able 
to treat these systems more accurately. QMC also has 
the property of scaling well with system size, although 
with a large prefactor, and has seen limited application 
to TM solids. As more computing power becomes avail- 
able, QMC calculations of solids will become routine, and 
this study offers some insight as to the accuracy that will 
be achieved. 

In this paper, we expand our treatment to the first 
five TM-0 molecules (Sc,Ti,V,Cr,and Mn), studying not 
only the binding energy, but also the bond length and the 
dipole moment. To obtain the dipole moment, we apply 
the relatively new Reptation Monte Carlo0 method for 



the first time to heavy elements. We find that for binding 
energy and bond lengths, QMC offers unmatched accu- 
racy, while the dipole moment is in less agreement with 
experiment. We investigate the effect of the fixed-node 
condition on the dipole moment and find that while there 
is a significant nodal error, it is not enough to reconcile 
the calculation with the experiments. 



II. METHOD 

A. Quantum Monte Carlo 

We use the Variational, Diffusion, and Reptation 
flavors of Quantum Monte Carlo (VMC,DMC, and 
RMC) in our calculations as implemented in the QWalk 
program^. In VMC, we start with a Slater determinant 
of one-particle orbitals, ^hf-, or a linear combination 
of Slater determinants. We then multiply ^ hf by the 
explicitly correlated inhomogencous Jastrow correlation 
factor e u . We write 



(1) 



where the lower case indices stand for electronic coor- 
dinates, and the upper case indices are ionic coordi- 
nates. The correlation factor is expanded in the Schmidt- 
Moskowitz formQ: 

u(r u , r 3l ,r l3 ) = ^ c e k l a k (r u ) + ^ C& fc (r 4j ) 
k m 

+ c kTm i r u )m [r 3 i ) + a k {r 3 / )a t (r tI ) ) b k {r i3 ) , 

kim 



where the a k and b k functions are written as 

1 - z(r/r cut ) 



ak{r) = 



1 + /3z(r/r cut )' 



(2) 



The polynomial z(x) = x 2 (6 — 8x + 3x 2 ) is chosen so 
the functions go smoothly to zero at r cut =7.5 bohr. Wc 
generate random samples in the 3-/Ve-dimensional space 
(denoted by R) according to the many-particle probabil- 
ity distribution *(R) 2 . The energy is then obtained by 
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averaging the local energy £x(R) = ^wrT ■ Following 
the variational theorem, we then optimize a combina- 
tion of energy and variance of the local energy, using the 
algorithm proposed by Umrigar and Filippi|9|. (3 and 
all the coefficients are variationally optimized. We use 
the VMC wave function as a trial function for Reptation 
Monte Carlo or Diffusion Monte Carlo. 

DMC and RMC are based on the so-called imaginary 
time Schrodinger equation 



cZW(R,- 



= (H-Bb)*(R,r), 



(3) 



which has a steady-state solution when '5 is an eigenvalue 
with value Eq and all non-steady-state solutions converge 
exponentially to the eigenstate <&o as r goes to infinity. 
Transforming to an integral equation, we have 

$o(Ri)=lim /"G(Ri,Ro,r)* r (Ro)dRo (4) 

where G is the Green's function of the imaginary time 
Schrodinger equation and \E , t(Ro) is the trial wave func- 
tion that we obtain from VMC. Solving for the exact G 
for large r is as difficult as solving for $o, so we choose 
some constant small value of r for which we know G ac- 
curately, and compound the operations (suppressing the 
r dependence of G) : 

$o(R)= lim G(R,R„)...G(Ri,R )*t(Ro). (5) 

n — >oc 

Each application of G is interpreted as a stochastic pro- 
cess, in the same way that the diffusion equation can be 
mapped onto Brownian particles and vice versa (in fact, 
for a free particle, the Hamiltonian is — ^V 2 and the sim- 
ulation is a diffusion process). 

DMC performs a simulation of random particles for 
large n. Skipping over some details that can be found 
in Ref [ijj , we eventually find that it obtains Pr^ (R) = 
$o(R)^t(R), which can be used to evaluate the ground- 
state energy as follows: 



(E Q ) = J rfR* T(R)(i > o( R)_lL_2, 



(6) 



since $o is an eigenstate of H and H can operate for- 
wards or backwards. Any operators that do not com- 
mute with the Hamilonian will have expectation values 
in error. We account for this error by using Reptation 
Monte Carlo, where the random walk is performed in the 
space of paths: s — [R , Ri, . . . , R n -i, R«]- We sample 
the path probability distribution 

n(a) = * T (Ro)G(Ro, Ri) • . . G(R«-i, R„)*t(R«) 

This can be interpreted in several different ways. If we 
examine the distribution at Ro, we can view the samples 
of Green's functions as acting on $t(R h ), and therefore 
Pr (R ) = , I't(Ro) ( I ) o(Ro)- This is the same distribu- 
tion as we obtain DMC as the path length goes to in- 
finity. Alternatively, since G is symmetric on exchange 



of the two R coordinates, the probability distribution of 
R n is the same. Finally, we can split the path in two, 
one projecting on \E , t(Ro), and the other projecting on 
* T (R„). We then have 

Pfl„ /2 (R n/2 ) = (G(R n/2 ,R n/2 _ 1 )...G(Ri,R 2 )*r(Ro)) 
x(G(R„ /2 ,R„ /2+1 ) . . .G(R n _i,R„)* T (R„)) 

= $o(Rn/2) 

for n — > 00, which allows us to obtain correct expecta- 
tion values of operators that do not commute with the 
Hamilonian. 

We use a Diffusion Monte Carlo algorithm very simi- 
lar to that described in Ref UXj- The Reptation Monte 
Carlo algorithm is from Ref [llj, except that we use the 
approximation to the Green's function as described in 
Ref [13 • 

Diffusion Monte Carlo and Reptation Monte Carlo 
both suffer from the sign problem for fermions, which 
forces us to make the fixed-node approximation, where 
the nodal surface of the exact wave function are assumed 
to be the same as the trial wave function. This approxi- 
mation typically results in recovering 90-95% of the corre- 
lation energy. This and the pseudopotential localization 
approximation^^ are the only uncontrolled approxima- 
tions in our calculations. All calculations will be done 
using these two approximations. 

We can control the fixed node approximation some- 
what by varying the orbitals in the trial wave function 
and minimizing the DMC energy. For transition metal 
oxides, this turns out to be important, since Hartree-Fock 
orbitals are rather inaccurate and biased towards more 
ionic picture of the state. We have previously optimized 
the mixing percentage in B3LYP in Ref |5j, and it turns 
out that B3LYP orbitals are almost optimal, so in these 
calculations we simply use the B3LYP orbitals. We will 
also investigate using multiple determinants to improve 
the fixed node error for a case study of TiO. 



B. Bayesian Fitting of Bond Lengths 

The main disadvantage of QMC methods is that ev- 
ery quantity has a statistical uncertainty which decreases 
only as the square root of the computer time. For quan- 
tities like bond lengths, researchers have historically cal- 
culated the energy at several bond lengths, then fitted 
a function to the points. Uncertainties have been cal- 
culated in many ways, but to our knowledge, none of 
them is exact and makes use of all the information avail- 
able including the statistical uncertainty. Here we offer a 
more systematic way of finding the minimum bond length 
along with its statistical error bar. 

According to Bayes' theorem, given a model M and a 
set of data D, the probability of the model given the set 
of data is: P(M\D) = P(D\M)P(M)/P(D). P{D) is 
an unimportant normalization constant, P(M) is called 
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the prior distribution, which we are free to set to re- 
flect the a priori probability distribution on the set of 
models. Without any good reason to believe otherwise, 
we generally set P(M) = 1, the maximum entropy/least 
knowledge condition. In the case of normally distributed 
data on a set of points {xx, X2, xn}, 

P{D\M) oc exp[- ^2(M(xi) - D( Xi )) 2 /2a 2 ( Xi )}, (8) 

i 

where u(x) is the statistical uncertainty of D(x). 

In the case of our bond lengths, we limited our space 
of models to M(x) — c\ + c 2 x + c 3 x 2 , for x close to 
the minimum bond length. This is equivalent to setting 
the prior distribution equal to one for all quadratic func- 
tions and to zero for non-quadratic functions. We then 
calculated several data points D{x) with statistical un- 
certainties a(x). The probability distribution function of 
the bond length b was then calculated by evaluating the 
integral 

= J S(- Cl /2c 2 - b)P(M\D)P(M)d Cl dc 2 dc 3 
P[ ' J P(M\D)P(M)dc 1 dc 2 dc 3 ' [ ' 

This integral is only three-dimensional, and as such could 
be calculated by a grid method, but we found it con- 
venient to calculate it by Monte Carlo, by sampling 
P(M\D)P(M) and binning the bond length. In all cases 
studied, p(b) was very accurately a Gaussian distribu- 
tion function, and so we report the bond lengths as an 
expected value plus a stochastic uncertainty, which fully 
characterizes the distribution. 



C. Computational Parameters 

For the oxygen atom, we used the pseudopotential 
from Lester[l3|, and for the transition metals, we used 
Ne-core soft potentials from LeeQ. To prepare the or- 
bitals for the QMC calculation, we used GAMESS[I3 
with a triple-zeta basis set. Both RMC and DMC calcu- 
lations were performed with r = 0.01 Hartree , which 
was converged within error bars, and for our RMC cal- 
culations, we chose N = 301, which corresponds to a 
3 Hartree -1 long projection length. We evaluate the 
dipole moment within RMC as the expectation value of 
A 4 = e (J2i er i) + ^nuclei, using the Hellmann-Feynman 
theorem. 



III. RESULTS AND DISCUSSION 
A. Energetics 

We begin with the importance of the one-particle or- 
bitals used in the trial function. These are not optimized 
within VMC and the Jastrow factor does not change the 
nodes, so we are forced to use the nodes of the Slater 
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FIG. 1: The energy gain in DMC from using B3LYP orbitals 
as a function of the metal monoxide. The line is a guide to 
the eye. 



determinant of orbitals from DFT or Hartree-Fock. For 
systems without strong electron correlation(for example, 
first and second row elements), it has been standard prac- 
tice to use DFT and Hartree-Fock orbitals interchange- 
ably, and the fixed-node energy is fairly insensitive. In 
TMO's, the correlation is much stronger and changes the 
structure of the one-particle orbitals. That is, the or- 
bitals that define the lowest energy nodal structure are 
significantly different from the Hartree-Fock orbitals. Di- 
rect optimization of the orbitals within QMC is desir- 
able, but very difficult for larger systems, so we took 
the approach of finding an optimal mean-field that pro- 
duces orbitals that minimize the fixed-node energy. In 
these systems, the hybrid functional B3LYP appears to 
be near-optimal. In Fig^ we report on the energy gain 
in DMC by using B3LYP orbitals. In our five molecules, 
there are roughly three levels of energy gain, correspond- 
ing to the type of bonding. ScO has only one d electron in 
a a state, TiO and VO are respectively ct 1 ^ 1 and a 1 5 2 1 
and CrO and MnO are crM 2 ^ 1 and cr 1 ^ 2 ^ 2 . Each new 
type of symmetry adds approximately 0.2 eV to the en- 
ergy gain in using B3LYP orbitals, with a slight decrease 
for the half-filled shell of MnO. This energy gain is a 
measure of how poor the independent-electron approxi- 
mation is for preparing the one-particle orbitals. Since 
there is almost no gain in the atomic systems by using 
B3LYP orbitals, the correct orbitals appear to be criti- 
cal for high accuracy in TMO materials, even more so as 
more d-symmetry electrons are present. 

We compare our QMC results to Density Functional 
Theory in the LDA, Coupled Cluster with singles, dou- 
bles, and perturbative triples (CCSD(T)), and a new hy- 
brid meta-GGA, TPSSh, which should be the most ac- 
curate semi-empirical DFT available [l^. Using accurate 
one-particle orbitals, DMC binding energies (Table HJ all 
fall within experimental uncertainties except for CrO and 
MnO, which both have 7r-type electronic configurations. 
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Method 


ScO 


TiO 


VO CrO MnO 


RMS 


LDA[1] 


9.09 


9.13 


8.48 6.26 6.51 


2.19 


CCSD(T)[2] 


6.71 


6.64 


6.13 4.20 3.43 


0.31 


TPSSh[l] 


7.11 


7.18 


6.44 4.45 4.62 


0.38 


DMC 


7.06(3) 


6.81(3) 


6.54(3) 3.98(2) 3.66(3) 


0.21 


Expfl7] 


7.01(12) 


6.92(10) 


6.44(20) 4.41(30) 3.83(8) 






TABLE I: Binding energies of the first five transition metal 
monoxides by different theoretical methods, along with RMS 
deviations from the experiment (all in eV). Statistical uncer- 
tainties in units of 10 -2 eV are shown in parentheses for 
Monte Carlo and experimental results. Zero point energy 
corrections are estimated to be much less than the size of 
the uncertainty in experiment. 



Method 


ScO 


TiO 


VO 


CrO 


MnO 


RMS 


LDA[1] 


1.644 


1.597 


1.564 


1.584 


1.602 


0.033 


CCSD(T)[2] 


1.680 


1.628 


1.602 


1.634 


1.66 


0.011 


TPSSh[l] 


1.659 


1.613 


1.582 


1.612 


1.628 


0.012 


DMC 


1.679 


1.612 


1.587 


1.617 


1.652 


0.008 


Expfl7] 


1.668 


1.623 


1.591 


1.621 


1.648 





TABLE II: 


Bond 


lengths 


in A. 


The statistical 


uncer- 



tainties for ScO,TiO,VO,CrO, and MnO are respectively 
0.002,0.003,0.003,0.004, and 0.004. 

The RMS deviations of DMC are around 50% smaller 
than TPSSh and CCSD(T) at 0.21 eV. This is still above 
the systematic error of 0.05 eV that would be required 
for 'chemical accuracy'; however, the uncertainties of the 
experiments are also above this threshold. 

TablelTTlshows the calculated versus experimental bond 
lengths for the selected methods. Here again, we see 
that DMC using a Slater determinant of B3LYP orbitals 
is quite accurate, with RMS deviation less than 0.01 
A, close to the size of our statistical error. It is again 
around 50% more accurate than the high-accuracy meth- 
ods CCSD(T) and TPSSh, and four times more accurate 
than the LDA on these systems, on average. 



B. Dipole Moment 

The dipole moment of these molecules has been noted 
as a difficult quantity to accurately calculate 0,0- In ap- 
proaches relying on basis functions, there appears to be 
a large sensitivity to quality of the basis set. It also ap- 
pears to be very sensitive to an accurate treatment of the 
correlation. The RMC method depends only very weakly 
on the basis set used to prepare the orbitals, and reaches 
the lowest energy of any variational method on these sys- 
tems, so one may hope that RMC agrees with experiment 
better than other ab initio methods. We find this not to 
be the case. As shown in Table IIIII we find serious dis- 
agreement with experiment in three of the four molecules 
with experiments available. Only ScO, the molecule with 
the weakest d-character, agrees well. The rest are uni- 



Method 


ScO 


TiO 


VO 


CrO 


MnO 


LDA[1] 


3.57 


3.23 


3.10 


3.41 




CCSD(T)f2] 


3.91 


3.52 


3.60 


3.89 


4.99 


TPSSh|l] 


3.48 


3.43 


3.58 


3.97 




RMC 


4.61(5) 


4.11(5) 


4.64(5) 


4.76(4) 


5.3(1) 


Exp [18] 


4.55 


3.34(1)[19] 


3.355 


3.88 





TABLE III: Dipole moments in Debye. The fixed-node RMC 
results have been obtained with a single deteriminant of 
B3LYP orbitals. See text for an analysis of the errors in- 
volved for the case of TiO. 

versally predicted to be much higher in QMC. 

The significant discrepancies from experiment are sur- 
prising given the excellent agreement that we obtained 
with energies. It also seems strange that the LDA is gen- 
erally quite close to the experiment, since we know that 
for energetics it performs relatively poorly. We analyze 
the errors present in the calculations as follows for the 
case of TiO, the simplest of the molecules with a large 
difference from experiment. 

• Pseudopotential error. We checked the pseu- 
dopotcntial in mean-field calculations, and it 
caused an overestimation of the dipole moment in 
TiO by 0.1 Debye. This is systematic for all five 
materials studied, with each having an overestima- 
tion of between 0.1 and 0.15. 

• Hellmann-Feynman theorem. The definition of 
dipole moment is = d j^) , where E is the electric 
field. We have used the Hellmann-Feynman the- 
orem to instead evaluate it as j i = e(^Yer;} + 
l^nuciei- As shown in Ref [2fJ, the Hellmann- 
Feynman theorem for calculating the dipole mo- 
ment does not exactly apply in fixed-node Quan- 
tum Monte Carlo, although the errors are generally 
considered small. To check this, we calculated the 
dipole moment usin g th e finite field approach and 
correlated sampling[2lj using nodes from B3LYP 
also at that field, and obtained the same result 
as the Hellmann-Feynman estimator within error 
bars. 

• Fixed-node error and localization. The only 
remaining approximation in our simulation is the 
combined fixed-node error and localization error. 
We investigate this by attempting to systematically 
improve the trial wave function. 

To go beyond the Slater-Jastrow form, we began to add 
more determinants into the trial function for the test 
case of TiO. We performed a configuration interaction 
calculation including single and double excitations start- 
ing from the B3LYP orbitals, kept the determinants with 
the highest weights, and reoptimized the weights in the 
presence of the Jastrow factor. If we did not reoptimize 
the weights, we found that the fixed-node energy actu- 
ally increases, suggesting that the correlation present in 
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the Jastrow factor is critical to an accurate description 
of these materials. The Jastrow factor also reduces the 
number of determinants necessary to describe the elec- 
tron correlation by several orders of magnitude, since the 
Jastrow factor contains most of the so-called dynamic 
correlation. Finally, we used an RMC simulation to find 
the dipole moment using the Hellmann-Feynman theo- 
rem. Fig |21 shows significant improvements in the energy 
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FIG. 2: The number of determinants versus the energy and 
dipole moment for TiO.The dipole moments are shifted down- 
wards by 0.1 Debye to correct for the pseudopotential error. 

on the order of .015 Hartrees, and the dipole moment ap- 
pears to oscillate around a value of approximately 3.9(1) 
Debye. Therefore, we estimate the error in dipole mo- 
ment to be approximately 0.2 Debye from the fixed-node 
approximation and 0.1 Debye from the pseudopotential 
approximation. With these corrections, the minimum 
value of the dipole moment for TiO is 3.6 Debye with over 
95% confidence. This is consistent with the CCSD(T) 
number, but still visibly larger than the value reported 
by experiment. 

IV. CONCLUSIONS 

We have found that for energetics, DMC using a single 
determinant trial function is remarkably accurate, per- 
haps suggesting that for these materials, it is sufficient. 



The bond lengths and binding energies are, on aver- 
age, 50% better than the best meta-GGA and CCSD(T). 
The Bayesian method for finding the minimum bond 
lengths mitigates the inconvenience of statistical uncer- 
tainty, while improving the performance by using all pos- 
sible information. 

The dipole moment appears to be more challenging, 
and requires a complicated treatment of the wave func- 
tion nodes to obtain a stable value with respect to 
changes in the nodes. We have obtained a converged 
value for TiO, however, and it is still somehwat higher 
than suggested by experiment. This is in line with the 
large values for the dipole moment obtained by CCSD(T) 
and B3LYP; agreement with CC method is particu- 
larly reassuring. In addition, there are to our knowl- 
edge only two experimental measurements of the dipole 
moment 1 191 ,22] from the same group, which report signif- 
icantly different moments(2.96 versus 3.34 Debye). This 
may indicate a sizeable uncertainty in the experiment 
as well. The apparent agreement of LDA is almost cer- 
tainly fortuitous, because LDA underestimates the bond 
length, which in turn causes the dipole moment to be too 
small. In fact, one would generally expect LDA to un- 
derestimate the dipole moment even at the correct bond 
length, since it tends to make the molecule not ionic 
enough. This may also indicate an inaccuracy in the 
experiment. One would expect that the corrections for 
the fixed-node approximation for VO, CrO, and MnO 
are similar or slightly greater than TiO, and therefore 
are around 0.2-0.4 Debye. It is well-known that dipole 
moment is particularly sensitive to contributions from 
single excitations which, however, have only very minor 
contributions to energy by coupling through doubles. It 
is therefore plausible that the optimization procedures 
which we have employed were not able to recognize the 
weak signal from singles and therefore the wavefunctions 
are still not perfect in this respect. The dipole moment 
remains an extraordinarily sensitive quantity that is a 
stringent test of theory and experiment. It may be inter- 
esting to see if there are other wave functions that can 
describe the nodal surface to sufficient accuracy while be- 
ing more compact than a large determinantal expansion. 
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